Critical packing in granular shear bands 
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In a realistic three-dimensional setup, we simulate the slow deformation of idealized granular 
media composed of spheres undergoing an axisymmetric triaxial shear test. We follow the self- 
organization of the spontaneous strain localization process leading to a shear band and demonstrate 
the existence of a critical packing density inside this failure zone. The asymptotic criticality arising 
from the dynamic equilibrium of dilation and compaction is found to be restricted to the shear band, 
while the density outside of it keeps the memory of the initial packing. The critical density of the 
shear band depends on friction (and grain geometry) and in the limit of infinite friction it defines a 
specific packing state, namely the dynamic random loose packing. 

PACS numbers: 45.70.Cc, 81.40.Jj 
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I. INTRODUCTION 

Packing density of different particulate systems is of 
main interest for scientific fields including, but not lim- 
ited to, suspensions, metallic glasses, molecular systems, 
and granular materials. In three dimensions, for identi- 
cal spheres, the face-centered cubic (FCC) packing is the 
maximum possible [1|. This fills the space with a vol- 
ume fraction of 7r/(3'\/2) ~ 0.74. Random arrangements 
have much lower densities [3| • Different experiments and 
computer simulations revealed that the largest obtainable 
volume fraction of a random packing of identical spheres 
is around 0.64. This is known as the random close pack- 
ing (RCP) limit. A mathematical definition of this limit 
can be given through the concept of maximally random 
jammed state [SJ, 14| . 

Reaching the RCP limit needs careful preparation (e.g. 
tapping and compression). If glass or marble beads are 
simply poured into a container the volume fraction is usu- 
ally only around 0.6. A random, loose packing (RLP) at 
its limit of mechanical stability obtained by immersing 
spheres in a fluid and letting them settle has a volume 
fraction of 0.555 [5|. The volume fraction of RLPs ob- 
tained with different methods (both experimental and nu- 
merical) show that this packing state is less well defined 
than the RCP limit. Attempts made in order to relate 
RLP to rigidity percolation [5| and to critical density at 
jamming of an assembly of (infinitely) rough spheres [a], 
are to be mentioned. 

Already in 1885, Reynolds noted that dense granular 
samples dilate during slow deformation [3] . On the other 
hand, it is well known that loose granular materials den- 
sify in such a process [^, Q • Under slow shear the strain is 
usually localized to narrow domains called shear bands. 
As it was first suggested by Casagrande |Tol , it is tempt- 
ing to assume that in these failure zones the system self- 
organizes its packing density to a critical value indepen- 
dent of the initial packing state of the material. 



While this hypothesis forms the basis of many contin- 
uum constitutive models of soil mechanics since decades 
jllj l , a general micromechanical theory of shear band for- 
mation and of the involved criticality is still missing. 
Progress, needed in order to deepen our understanding 
of the critical state in shear bands, can be expected from 
the remarkable development of experimental techniques 
(including Computer Tomography [12l . [l3| and measure- 
ments in microgravity |14l . Il5|) and of simulation tech- 
niques which become increasingly efficient as computa- 
tional power grows |16j . 

Shearing of granular materials has been investigated 
in many different geometries and specially designed lab- 
oratory tests (for recent results see [l2, llJl)- Such ex- 
perimental studies revealed complex localization patterns 
and presented evidence for the existence of a critical par- 
ticle density inside the shear bands. The importance of 
computer simulations is enhanced by the fact that they 
make possible studies which are difficult to control in ex- 
periments (e.g. friction dependence) and they facilitate 
the measurement of hardly accessible quantities (e.g. vol- 
ume fraction inside the shear bands). 

The critical density, in numerical studies, is often stud- 
ied only in special conditions when shearing extends to 
the whole volume of the samples. This allowed for dis- 
cussing the criticality based simply on global behavior 
(e.g. dilatancy). Without reference to shear bands, many 
qualitative effects were already pointed out in both math- 
ematical models [13, ^M ^-i^d simulations [13, ll^ • How- 
ever, such studies neglected the involved localization phe- 
nomena inevitable in real situations and disregarded the 
self-organizing manner in which the packing state of the 
shear bands is usually formed. 

A principal parameter which controls the dynamic 
equilibrium between dilation and compaction in fully de- 
veloped shear bands is the friction between the grains. In- 
tuitively, a system of frictionless grains can be sheared at 
a large packing density (close to the RCP limit) because 
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FIG. 1: (Color online) a) Grains placed between two hori- 
zontal platens and surrounded by an elastic membrane were 
subjected to a vertical load and a lateral confining pressure, 
b) The membrane was modeled with overlapping spheres [20| . 



the grains (under slow shear) can easily rearrange in com- 
pact configurations. At large friction the rearrangement 
of the grains is hindered by friction, consequently the 
packing density of the shear bands is expected to define 
a low density state close to the RLP limit. 

The aim of this Paper is to study the emergence of a 
critical packing state in sheared granular media and to 
present its relation to shear bands as well as its depen- 
dence on friction. 



II. SIMULATIONS 

We investigate numerically an axisymmetric triaxial 
shear test (see Fig. [T|). This consists of the slow com- 
pression of a cylindrical sample enclosed between two 
end platens. The sample is surrounded by an elastic 
membrane on which an external confining pressure is ap- 
plied. The end platens are pressed against each other in 
a strain controlled way. The upper platen is allowed to 
tilt. In certain conditions a planar shear band is formed 
[12l | . Using different initial packing densities and friction 
properties of the grains as well as identifying the grains 
in the failure zones makes it possible to study the critical 
packing density inside the shear bands. 

The simulations, which are going to be presented here, 
are based on a standard three-dimensional Distinct El- 
ement Method (DEM) [2l|. We implemented the Hertz 
contact model [22l | with appropriate damping combined 
with a frictional spring-dashpot model [16| . The normal 
Fn and the tangential Ft forces are calculated as 
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where K„ = 10^N/m^'% Kt = 10'^N/m, 7„ = lNs/m' 
and 7t = 1 N s/m are normal and tangential stiffness and 
damping coefficients, (5„ and 8t are normal and tangen- 
tial displacements, and v^ and Vj are the normal and 
tangential relative velocities. 



TABLE I: Volume fraction r\Q of samples prepared with dif- 
ferent coefficients of friction ^o • With each ^o we prepared 2 
samples having the same 770 within 0.2% relative error. 



The numerical values are chosen to realize the hardest 
material that we could safely simulate with the minimal 
damping preserving the numerical stability of the calcula- 
tions. With the above stiffness and damping coefficients, 
the inverse of the average eigenfrequency of contacts, in 
both normal and tangential direction, is more than one 
order of magnitude larger than the used integration time 
step At = 10~^ s. This assured that the noise level in- 
duced by numerical errors is kept low. 

With relatively small samples (made up of 27000 par- 
ticles) but in a realistic geometry we have succeeded to 
reproduce shear band morphologies [20, [23| known from 
experiments [IJ, ll3|- Iii order to study the criticality 
of these shear bands, we prepared homogeneous initial 
configurations of different volume fractions using the de- 
position method described in |20| . 

We used a particle distribution similar to those en- 
countered in experimental studies of idealized granular 
materials. Our particles are spherical, they have equal 
mass density (2.5 • 10'^ kg/m ), equal friction coefficient, 
and their diameters are set according to a narrow Gaus- 
sian distribution with mean d = 0.9 mm and standard 
deviation of 2.77%. The prepared cylindrical samples, 
having diameter D = 23.3 d, consisted of 20000 to 27000 
spherical grains as required by a prescribed packing den- 
sity and the H k,2.2 D geometrical constraint, where H 
is the height of the samples. 

Initially the particles were placed randomly in a tall 
cylinder (about 3 times taller than H). They were given 
small downwards velocities in such a way that they all col- 
lided approximately at the same time. The upper platen 
was pressed on top of the packing to hold it together. 
This method provides an efficient way to produce a ho- 
mogeneous random packing. The volume fraction of the 
prepared samples could be controlled in the full RLP to 
RCP range (see Tab. |T]) by varying the coefficient of fric- 
tion /iQ which was applied during this phase. 

After preparation, the friction coefficient of the parti- 
cles was set to a new value /x independent of ^o- During 
the simulations, similarly to Q, we compressed the sam- 
ples vertically at zero gravity and 0.5 kPa confining pres- 
sure. The bottom platen was fixed. The upper platen 
moved downward with a constant velocity, inducing an 
axial strain rate of 20 mm/s. During compression the up- 
per platen could freely tilt along any horizontal axis with 
rotational inertia / = 10^^ kg m^ . 

The lateral membrane surrounding the sample was 
modeled with approximately 15000 identical, overlap- 
ping, non-rotating, frictional spheres connected with elas- 
tic springs. The stiffness of the springs was set to 



Ks = 0.5 N/m. This prevented the partieles from escaping 
by passing through the membrane. The membrane parti- 
cles were initially arranged in a triangular lattice (Fig. [1] 
(b)). The confining pressure was applied on the triangu- 
lar facets formed by the neighbo ring "m embrane nodes" 
as described in (see also 0, [iSllil ) . 

It is worth mentioning that Cui and O 'Sullivan j26| 
have introduced a technique which speeds up calculations 
by computing only a section of the cylindrical sample. 
This allows for larger samples but requires that the sym- 
metry of the system is kept during compression, and thus 
eliminates the possibility of symmetry breaking strain lo- 
calization, which arises spontaneously [20| if tilting of the 
upper platen is not suppressed (see Fig. [T](a)). 

We have executed several simulation runs. The grain- 
platen and grain-membrane contacts were calculated sim- 
ilarly to grain-grain contacts including the friction prop- 
erties. The samples of different initial volume fractions 
(see Tab. |I| were first compressed using the same coef- 
ficient of friction /i = 0.5. Later, we compressed the 
densest samples (ryo = 0.641) with 10 different friction 
coefficients ^ G {0, 0.1, . . .0.9}. For each set of parame- 
ters, two simulation runs were executed using specimens 
prepared with different random seeds. 

During compression, we measured locally the shear in- 
tensity S and the volume fraction 77. The regular trian- 
gulation of the spherical grains [23, [2^ was used to de- 
fine these quantities. The local volume fraction is given 
by the ratio of the volumes of a grain and its regular 
Voronoi cell. The local shear intensity is calculated from 
the macroscopic strain tensor derived from particle dis- 
placements [20, [23] . Using the eigenvalues e^ of this ten- 
sor, we defined the local shear intensity as 



S = max 

k 



£k 



It.'^ 



(3) 



To overcome fluctuations due to random packing and re- 
arrangements, we calculated spatial averages up to 3rd 
order neighbors along the regular triangulation. 
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III. IDENTIFICATION OF HIGH SHEAR 
INTENSITY REGIONS 

Strain localization in dense and loose samples shows 
substantial differences [i2| • In dense samples shear bands 
are usually formed after a short plastic deformation and 
inside them the local packing density is lower than in the 
bulk (i.e. the regions outside of the shear bands). Since 
denser parts arc more stable, the position of the shear 
bands remains unchanged for the whole duration of a 
shear test. Contrary, in loose samples the shear bands 
have a slightly higher packing density than the bulk and 
hence the position of the shear bands is likely to change 
and to move around the whole sample. This leads to more 
or less homogeneous samples with local packing densities 
close to the packing density of the shear bands. 



FIG. 2: (Color online) Example of shear intensity and vol- 
ume fraction histogram maps (rjo — 0.641, jj, — 0.5). The 
shear intensity is measured in arbitrary units. On (a, b) the 
occurrences (scaled to [0, 1]) are encoded with the color scale 
shown at the top. On (a) the white curve marks the shear 
intensity threshold Sth ■ On (b) it marks the shear band vol- 
ume fraction rjsB- The dotted and dashed vertical lines on 
(a, b) at 4.4% and 16% axial strain mark the position of the 
histograms (c, d) and (e, f), respectively. 



The general algorithmic identification of failure zones 
based on geometric methods is difficult, especially, re- 
garding the identification of the non-persistent shear 
bands of the loose samples. Nevertheless, based on the 
calculated local shear intensity, individual grains can be 
categorized to be part of the failure zones or the bulk 



providing that a good enough threshold separating these 
two classes can be found. A histogram technique seems 
to be a perfect candidate for this. 

Fig. [5] presents histograms of the local shear intensity 
and local volume fraction. It shows the main aspects 
of the strain localization process observed in one of our 
simulations executed with a sample having rjQ = 0.641 
and fi = 0.5. At the beginning of the test, up to approxi- 
mately 6% axial strain, in shear intensity histograms (see 
Fig. [5] (a, c)), we find a single peak at medium S. This 
means that almost all particles rearrange simultaneously 
and consequently the sample experiences a more or less 
plastic deformation. This is underlined by the fact that 
the local volume fraction has just one strong peak (Fig. [2] 
(d)), indicating that the sample is still homogeneous. 

At higher (> 6%) axial strain a shear band is formed. 
This is localized to a planar failure zone of width of ap- 
proximately 10 particle diameters and is characterized by 
much higher S than the bulk. In the bulk the shear in- 
tensity fluctuations are small, while these fluctuations are 
large in the shear band. Consequently, in shear intensity 
histograms (Fig. [5] (e)), we find a narrow peak at low 5, 
corresponding to the bulk, and a wide peak at high S, 
corresponding to the shear band. The volume fraction 
histogram does also become more structured showing ev- 
idence of a non-homogeneous material. The narrow peak 
at low volume fraction corresponds to the shear band 
while the bulk produces a much wider distribution at a 
higher volume fraction (Fig. [2] (f)). 

Motivated by this separation, we computed a shear 
intensity threshold Sth, which could be used to define 
two classes of shear intensity values (low and high) and 
to classify the grains accordingly into shear band and 
bulk. For this we used Otsu's threshold selection method 
[30[ described in Appendix [X] This histogram technique 
minimizes the within-class variance and maximizes the 
separation of classes, and thus gives an ideal solution to 
our problem. We have also tested another threshold se- 
lection method modeling the histograms with the sum of 
two Gaussian functions, however, numerically this proved 
to be less stable and less reliable. 

The condition S > Sth, made it possible to identify 
the grains in high shear intensity regions - where shear 
bands emerge - and thus the average volume fraction tjsb 
of these regions could be calculated. The local packing 
density in shear bands is found to have small fluctuations 
and to give a peak in the volume fraction histograms. 
This coincides with tjsb (see Fig. [2] (b, f)), giving a self 
validation of the method. 

Let us note that even if Fig. [5] (e) does not suggest a 
clean separation of the grains into those within and those 
outside the shear bands ~ i.e. there is no shear intensity 
gap between the two regimes - this is not crucial for rjsB ■ 
Adding an artificial random noise of 10% to Sth does 
influence the resulting tjsb only within 0.5%. 

We also mention here, that before strain localization 
takes place the shear intensity histograms have only one 
peak (see Fig. [5] (c)). In this case, the threshold given 
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FIG. 3: (Color online) Volume fraction measured in high 
shear intensity regions (a) and globally (b) as function of the 
axial strain. The different lines correspond to different ini- 
tial packing densities rjo (see Tab. |T]) decreasing from top to 
bottom. The two panels use the same notations. 



by Otsu's method, which falls on the middle of the peak, 
is not physically relevant. However, as the sample is ho- 
mogeneous, the selected samples in high shear intensity 
regions still give the volume fraction which is close to the 
average volume fraction of the whole sample. This can 
be verified on Fig. [2] (d). 



IV. RESULTS 

As expected [13, [iJI , we found that due to strain lo- 
calization at the end of the shear tests the global volume 
fraction of the samples is not equal to the packing den- 
sity of the high shear intensity regions and thus global 
measurements cannot be used to characterize the prop- 
erties of failure zones. The behavior of both dense and 
loose samples demonstrates that in the shear bands, the 
initial packing conditions are canceled and a critical vol- 
ume fraction ry^ is reached in a self-organizing manner 
independently of the initial density of the tested granu- 
lar specimens (Fig. |3](a)). 

The criticality is found to be restricted to the shear 
bands. The global volume fraction rjg calculated from 
the total volume of the samples does not converge to 
rjc (Fig. [3] (b)). This behavior is expected to be more 
pronounced on larger systems. The dense samples are 
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FIG. 4: (Color online) Critical packing of shear bands as 
function of friction measured at 20% axial strain. The fitted 
curve shows t?c(m) = '??° ~ iv'H" ~ Vc) exp(— /i//x°), where rj° — 
0.637, ??r = 0.578, and fi° = 0.23. 



characterized by rjg > ijc- This demonstrates that the 
dilatancy [3 is concentrated to the shear bands. Con- 
trary, for loose samples rjg < rj^, however, the specimen 
is only slightly looser outside the shear bands. This gives 
a direct proof of shear induced compaction [9| . 

As we could see, at given friction, the critical pack- 
ing of shear bands is independent of the initial den- 
sity. It is, however, a further question whether it de- 
pends on friction. We studied quantitatively this effect 
for samples with initial density 770 = 0.641. At /i = the 
high shear intensity regions have a large volume fraction 
?7° = ?7c(0) = 0.637 ±0.002, which is only shghtly smaller 
than 770, showing that frictionless granular systems can 
be sheared at densities very close to the RCP limit. 

In frictional systems, the volume fraction of the fully 
developed shear bands is substantially lower than the ini- 
tial volume fraction (see Fig. [4]). With increasing /i, this 
decreases and converges to a limit which we estimate to 
77^ = lim^_^oo VcifJ-) = 0.578 ± 0.003 based on the expo- 
nential extrapolation of our data. 

This limit volume fraction depends only on geometry 
factors such as shape and size distribution of the grains 
and is characteristic to the dynamic equilibrium between 
dilation and compaction developed in a self-organized 
manner through strain localization. Based on the value 
of rj^ and this latter aspect, the corresponding asymp- 
totic state - which we refer to as the dynamic random 
loose packing (DRLP) ~ should be distinguished from the 
static RLP limit. 



V. CONCLUSIONS AND DISCUSSION 

We have presented distinct element simulations of ax- 
isymmetric triaxial shear tests at zero gravity and low 
confining pressure. Due to spontaneous strain localiza- 
tion, shear bands were formed. Using a histogram tech- 
nique, we identified the grains in high shear intensity re- 
gions, which at large axial strains coincide with the shear 



bands. We measured the packing density rjsB inside these 
failure zones and we found that in fully developed shear 
bands rjsB approaches a critical value 77c independent of 
the initial density of the samples. This in agreement 
with Casagrande's [lO| observation made for sandy soils 
seven decades before and also with recent experiments 
[13, [H, M, [II and numerical studies [13, [H, [ll . 

Rothenburg and Kruyt [l8| obtained similar results in 
two-dimensional simulations of biaxial shear tests. They 
have presented a theory of the average coordination num- 
ber of sheared granular media and derived a law for its 
evolution during slow deformations. Analyzing the re- 
lationship between volume fraction and average coordi- 
nation number, they conclude that a proper characteri- 
zation of granular media undergoing shear deformation 
should be based on packing density. 

We have shown that the criticality is restricted to the 
shear bands and global measurements (such as dilatancy) 
are unsuitable for the investigation of the properties of 
sheared granular materials in realistic situations, where 
strain localization is inevitable. To our knowledge, it is 
the first time that the critical packing density of shear 
bands was evidenced based on simulations of a realistic 
three-dimensional setup and spontaneous strain localiza- 
tion revealing the self-organizing manner in which the 
packing state of the shear bands is developed. 

We have further shown that 77c depends on the coeffi- 
cient of friction /z and in the limit /i — > cx) it converges to 
a value 77^ , which we have calculated within the accuracy 
of our simulations. The found 7]"^ defines a low density 
dynamic random loose packing (DRLP) state, which is 
characteristic to the dynamic equilibrium between dila- 
tion and compaction in the shear bands and depends only 
on the geometry of the grains. Based on the underlying 
mechanism, we argue that the asymptotic packing state 
of shear bands differs from the static RLP limit. 

This result should be also compared with the findings 
presented recently by Zhang and Makse [y| regarding the 
critical density of granular materials at jamming tran- 
sition. The importance of these findings lie in the fact 
that jamming is a basic concept which through a uni- 
fying phase diagram [3l|, [33] connects granular matters 
with a variety of other systems including dense particu- 
late suspensions and effects such as diverging viscosity at 
a maximum packing fraction [33, [3J| . 

In quasistatic limit, Zhang and Makse [6| reported a 
monotonous decrease of the critical packing density as 
function of friction. For low friction, the density of the 
shear bands found in our simulations is lower than at 
jamming found by Zhang and Makse [6|, while at high 
friction the situation is reversed. This indicates a natural 
separation of low and high density regions with possibly 
different mechanisms of dissolving the jammed state. 

As a final remark let us note that our results are de- 
rived for idealized granular materials composed of spheres 
having a narrow size distribution. It is well known that 
for non-spherical grains and wide size distributions the 
packing efficiency increases [35| , which should be also re- 



fleeted in the paeking density of the shear bands. This 
could be the reason why experimental results on sand re- 
veal smaller volume fractions in shear bands [12, 13] than 
the values found in our simulations. 



The within- class variance 



cr^it) = qiit)af{t) + q2{t)al{t) 



(A7) 
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APPENDIX A: OTSU'S THRESHOLD 
SELECTION METHOD 

Otsu's method [30| is a histogram technique known 
from Digital Image Processing, where it is typically used 
to transform grayscale images into two component (black 
and white) images. 

Let us consider a normalized histogram P{i), i.e. a 
histogram with the property 



E ^(*) = 1' 



(Al) 



where i is the bin index. The mean fj, and the variance 
a^ can be calculated as 



i 

a^ = Y.{^-^,fP{^). 
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^l (t) - qi {t)q2 [t) {fii (t) ~ 112 it)) ' ( A8) 

is a measure of the separation of classes. It is easy to show 
(T^. Otsu [3^ proposed to calculate 



that cr^ 



topt by either minimizing cr^ (t) 
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an optimal threshold t 
or maximizing ag{t). 

Maximizing cr^(i) is easier. It can be seen that 
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and thus 
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For each candidate threshold t, qi{t) and /ii(i) can be 
calculated with the recursive formula 
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Let us further consider a candidate threshold t and 
split the histogram in two parts Ti{t) = {i \ i < t} and 
X2{t) = {i\i>t}. With k G {1, 2} and 



qk 



.(t)= E ^W' 



(A4) 



iGJfc(t) 



the mean /^^(i) and variance cr^(t) of the two parts are 
defined by the equations 

qk{t)Mt) = E *^W' (^5) 

qk{t)al{t) = E (*-Mfc(t))'P(*). (A6) 
ieikW 



where gi(0) = P(0) and /ii(0) = 0. 

Both qi{t) and /ii(i) are increasing monotonously with 
t, consequently the maximum of (j^(i) is well-defined, ex- 
cept for degenerated cases which must be handled sepa- 
rately. The optimal threshold topt is given by the smallest 
candidate threshold s which satisfies the equation 
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Because ayy{topt) + o'B{topt) = cr^, the method both min- 
imizes the within-class variance and maximizes the sep- 
aration of classes. 
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